相関係数とは,数値データ同士の関連性を探る指標です.相関係数の絶対値が0に近いと2つの変数同士には線形関係がないことを示します.
ちなみに,この「相関の強さ」について分野によって評価が異なります.例えば,社会科学研究では高い相関が認められることは少ないです.今回の基準で中程度の相関や低い相関で議論をすることもあります.
この辺は分野によって異なりますので,ご承知おきください.
| dataset | 平均値 | 標準偏差 | 標本数 | 標準誤差 |
|---|---|---|---|---|
| away | 54.27 | 16.77 | 142 | 1.407 |
| bullseye | 54.27 | 16.77 | 142 | 1.407 |
| circle | 54.27 | 16.76 | 142 | 1.406 |
| dino | 54.26 | 16.77 | 142 | 1.407 |
| dots | 54.26 | 16.77 | 142 | 1.407 |
| h_lines | 54.26 | 16.77 | 142 | 1.407 |
| high_lines | 54.27 | 16.77 | 142 | 1.407 |
| slant_down | 54.27 | 16.77 | 142 | 1.407 |
| slant_up | 54.27 | 16.77 | 142 | 1.407 |
| star | 54.27 | 16.77 | 142 | 1.407 |
| v_lines | 54.27 | 16.77 | 142 | 1.407 |
| wide_lines | 54.27 | 16.77 | 142 | 1.407 |
| x_shape | 54.26 | 16.77 | 142 | 1.407 |
datasaurus<-datasaurus_dozen %>%
ggplot(aes(x=x, y=y, colour=dataset))+
geom_point()+
theme_void()+
theme(legend.position = "none")+
facet_wrap(~dataset, ncol=3)
datasaurus
x <- rnorm(100, 50, 10) y <- rnorm(100, 50, 10) z <- (x+y) / 2
plot(x, y)
plot(x, z)
\[ s_{xy}= \Sigma((xの偏差)*(yの偏差))/(データの個数-1)\]
x_hensa <- x-mean(x) y_hensa <- y-mean(y) goukeixy <- sum(x_hensa * y_hensa) kyobunsanxy <- goukeixy/(length(x)-1) kyobunsanxy
## [1] -7.868727
cov(x, y)
## [1] -7.868727
cov(x, z)
## [1] 43.82204
\[r = \frac{s_{xy}}{s_{x}s_{y}}\]
soukanxy <- kyobunsanxy/(sd(x)*sd(y)) soukanxy
## [1] -0.07328934
cor(x, y)
## [1] -0.07328934
cor(x, z)
## [1] 0.6333938
一般線形モデルとは,統計学の中でも,以下の数式(モデル式)を元に考えていくモデルです.
\[Y_i=\beta_1 X_1 +\beta_2 X_2+ \beta_3 X_3 + .... \alpha+\epsilon_i \]
さて,何か複雑そうなモデル式が出てきてしまいましたが,恐れることはありません.少し,簡単な形にしてあげましょう.そうすると,こんな感じに書くことが出来ます.
\[Y_i=\beta_1 X_1 + \alpha+\epsilon_i \]
このモデル式,何だか見覚えのあるグラフとそっくりだと思います.中学校の時に“一次関数”というのを教わったのを覚えていますでしょうか?一次関数ではこんな数式を使いました.
\[Y=\beta X + \alpha\]
この数式を元に,グラフを書く,ということもやったかと思います.この時,\(\beta\) を傾き,\(\alpha\) を切片という呼び方をしていました.ちなみに,この数式で直線のグラフを書く時には,Xに0を代入した時のポイント(0, \(\alpha\))とXに1を代入したときのポイント(1, \(\beta + \alpha\))を結ぶ直線を引いてあげれば,グラフを作成することができます.
一般線形モデルの一番理解しやすい最初の考え方は,「実際に観察されたデータを元にして,一次関数のような直線を引いてあげよう!」という発想です.ただし,一次関数とちょっと違うのは「全ての点を通らなくてよい」ということです.
一次関数の場合はその直線上にある全ての点を通ることが前提となっていました.しかし,実際には直線であるので,直線上の2点を通れば,全てその条件を満たす直線を引くことが出来ます.
しかし,一般線形モデルの場合は常に全ての点を通るとは限りません.ベストは全ての点を通ることではありますが,実際にはデータには「誤差」というものが存在します.これは本来得られるべき結果と実際に得られた結果にずれがあることを示しています.
この誤差には大きく分けて次の3種類あります.
測定誤差:実際に何かを計測する時に生じる誤差.中でも以下の2種類がある.
さて,少し本題に戻りましょう.ちょっと一般線形モデルのモデル式を考えたいと思います.
\[Y_i=\beta_1 X_1 + \alpha+\epsilon_i \]
改めて,このモデル式を説明したいと思います.ここで,“\(Y_i\)”のことを“応答変数”,“\(X_1\)”のことを“説明変数”と呼びましょう.
文字についている“\(_i\)”は各データによって異なる!という区別をするためについています.ちなみに,“\(Y_i\)”は他にも,被説明変数と呼ばれたりします.
また,"\(\beta_1\)" は係数,“\(\alpha\)” は切片と呼ばれます.そして,“\(\epsilon_i\)”が一番問題となる誤差です.この誤差は予測されたモデル式である“\(Y_i=\beta_1 X_1 + \alpha\)”からどれだけそのデータの値が離れているかを示しています.
と,言ってもなかなか理解し難いと思うので,一つ試しにやってみましょう.ここでは,「回帰分析」という方法と「t検定」という方法についてお話をしたいと思います.
| 検定名 | 応答変数 | 説明変数 |
|---|---|---|
| 回帰分析 | 数値データ | 数値データ(順序データ) |
| t検定 | 数値データ | 因子データ(ダミー変数,1, 0) |
回帰分析とは,応答変数が数値データであり,説明変数も数値データである場合に用いる方法です.例えば,「身長」と「体重」の間の相関関係について分析をする際にも用います.ここでは,今まで授業で使ってきた「主観的幸福度」と「生活満足度」の間に相関関係があるかどうか,以下の順番に沿って考えてみましょう.
この関係はモデル式で表すと,このような形になります.
\[(主観的幸福度)=\beta_1 (生活満足度) + \alpha+\epsilon_i \]
この時,切片であるは生活満足度が0であった時に対応する主観的幸福度を示しています.
何はともあれ,統計分析をするときには仮説を立ててあげる必要があります.仮説を立てるときには,「帰無仮説」と「対立仮説」の2つを考える必要があります. 対立仮説は「イイタイコト」,帰無仮説は「イイタイコトではないこと」でした.
ここで主観的幸福度と生活満足度の関係ですので,以下のように設定できます.
特に,以下では応答変数を主観的幸福度,説明変数を生活満足度とします.
はじめに,分析対象となるデータを読み込んでおきましょう. * もちろん,既に読み込んである場合は飛ばしてもらって構いません.
散布図のプロットは他の機能から持ってきてもよいのですが,今回はRStudio上でクリックだけで入れられる方法を紹介します.
その上で,コードを貼り付けて出力することにしましょう.
library(ggplot2) ggplot(exdataset) + aes(x = SUB_HAP, y = SUB_SAT) + geom_point(size = 1L, colour = "#0c4c8a") + geom_smooth(span = 1L) + theme_minimal()
どうもグラフを見ている限りだと,この2変数間には正の相関関係,すなわち「生活満足度が高ければ高いほど,主観的幸福度が高くなる」という傾向にはありそうです.
ただし,今はグラフを見ているだけなので,果たしてこの傾向が本当にあるのかどうかがわかりません.今度はこの傾向が科学的に認められるのかどうかを考えてみましょう.
さて,今度はRで分析してみましょう.ここでは,2行ほどのコードを書いてもらいます.
hapsat_model<-lm(SUB_HAP~SUB_SAT, data = exdataset) summary(hapsat_model)
## ## Call: ## lm(formula = SUB_HAP ~ SUB_SAT, data = exdataset) ## ## Residuals: ## Min 1Q Median 3Q Max ## -7.8918 -0.6503 -0.0814 0.7289 6.4015 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.59853 0.10176 15.71 <2e-16 *** ## SUB_SAT 0.81036 0.01711 47.37 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 1.285 on 961 degrees of freedom ## Multiple R-squared: 0.7002, Adjusted R-squared: 0.6999 ## F-statistic: 2244 on 1 and 961 DF, p-value: < 2.2e-16
## Call: ## lm(formula = SUB_HAP ~ SUB_SAT, data = exdataset)
この行では,分析したモデル式について示しています.簡単に言うと,「生活満足度によって,主観的幸福度は説明できるかどうか試してます...」ということを示しています.
## Residuals: ## Min 1Q Median 3Q Max ## -7.8918 -0.6503 -0.0814 0.7289 6.4015
ここでは,モデル式からのズレ(\(\epsilon_i\))である誤差がどの程度あるのかを示しています.ここでは誤差の最小値,第1四分位点,中央値,第3四分位点,最大値を示しています.一般線形モデルではこの誤差が正規分布になっていることを仮定しています.
## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.59853 0.10176 15.71 <2e-16 *** ## SUB_SAT 0.81036 0.01711 47.37 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
したがって,この結果はモデル式で書くと,以下のように示すことが出来ます.
\[(主観的幸福度)=0.81036 ×(生活満足度) + 1.59853+\epsilon_i \]
このモデル式は生活満足度が1あがると,主観的幸福度が0.8106ポイント増加すること,そして生活満足度が0である人の主観的幸福度は1.59853であることが推定されています.
ここに出てくるt valueはt値を,Pr(>|t|)はp値を示しています.そして,最後のsign.if. codesでは,どのような基準で*をつけているかを説明しています.この場合,p値が1-0.1の場合は無印,0.1-0.05の場合は“.”,0.05-0.01の場合は“*”,0.01-0.001の場合は“**”,0.001-0の場合は“***”,としてつけている,ということが示されています.
統計学の基本的な考え方ではp値が0.05以下,すなわち5%以下である場合には対立仮説を採択することがお約束となっています...が,単純に5%以下であることによって対立仮説を採択することがあってはいけません.
それは以下の理由によります.
## Multiple R-squared: 0.7002, Adjusted R-squared: 0.6999 ## F-statistic: 2244 on 1 and 961 DF, p-value: < 2.2e-16
続いて,確認したいのはこの2行です.“Multiple R-squared”はR2乗(あーるにじょう)値を示しています.ただし,このR2値は決定係数と呼ばれており,回帰式の当てはまり具合を示しています.寄与率とも呼ばれて,この値が1に近ければ近いほどよく説明できているモデル式であると言われます.ただし,R2乗値はこのモデルに組み込まれる説明変数が増えれば増えるほど,より良くなっていきます.そうするといくらでも興味のない変数を入れて重回帰分析(後日説明します)....と,なると決して意味があるモデル式になるとは言えません.
そこで,たくさん変数を入れたことに対するペナルティを加えたのが“Adjusted R-squared”,調整済みR2乗値と呼ばれるものです.こちらを報告してあげると良いかと思います.
最後の“F-statistic”はF検定と呼ばれるものの結果です.2つの群の「標準偏差」が等しいかどうか,を示しているものであり,「等分散性の分析」に用いられているものです.この結果は,主観的幸福度と生活満足度では分散,すなわちばらつき方が異なっている,ということを示しています.
結果の表記例.
生活満足度1が改善すると,主観的幸福度が0.81改善することが,0.1%水準で示された. (一緒に表を見せると良い.)
生活満足度1が改善すると,主観的幸福度が0.81改善することが示された(t(961)=47.37, p=.001).
\[(主観的幸福度)=0.81036 (t=47.37) ×(生活満足度) + 1.59853 + \epsilon_i \]
library(huxtable) huxreg(hapsat_model)
| (1) | |
|---|---|
| (Intercept) | 1.599 *** |
| (0.102) | |
| SUB_SAT | 0.810 *** |
| (0.017) | |
| N | 963 |
| R2 | 0.700 |
| logLik | -1607.061 |
| AIC | 3220.121 |
| *** p < 0.001; ** p < 0.01; * p < 0.05. | |
library(stargazer)
stargazer(hapsat_model, type = "html", align=TRUE,
title = "分析結果", out = "hapsatmodel.xls")
## ## <table style="text-align:center"><caption><strong>分析結果</strong></caption> ## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td><em>Dependent variable:</em></td></tr> ## <tr><td></td><td colspan="1" style="border-bottom: 1px solid black"></td></tr> ## <tr><td style="text-align:left"></td><td>SUB_HAP</td></tr> ## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">SUB_SAT</td><td>0.810<sup>***</sup></td></tr> ## <tr><td style="text-align:left"></td><td>(0.017)</td></tr> ## <tr><td style="text-align:left"></td><td></td></tr> ## <tr><td style="text-align:left">Constant</td><td>1.599<sup>***</sup></td></tr> ## <tr><td style="text-align:left"></td><td>(0.102)</td></tr> ## <tr><td style="text-align:left"></td><td></td></tr> ## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>963</td></tr> ## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.700</td></tr> ## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.700</td></tr> ## <tr><td style="text-align:left">Residual Std. Error</td><td>1.285 (df = 961)</td></tr> ## <tr><td style="text-align:left">F Statistic</td><td>2,244.149<sup>***</sup> (df = 1; 961)</td></tr> ## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr> ## </table>
\[t値=(期待値)-(平均値)/(標準誤差)\]
t値はこんな数式から算出されます.
標準誤差は(標準偏差)/(データ数の平方根)によって計算できることを思い出しておいて下さい. t値は分子が大きければ,平均値との差が大きいことを示しており,分母が大きければ,標準偏差(分散)が小さく,データ数が十分にあることを示しています.このt値が大きければ大きいほど,帰無仮説を棄却して対立仮説を採択できることを示しています.
一方,p値は帰無仮説が成立していることを前提として,0.05,すなわち5%未満であれば,帰無仮説を棄却するための基準となります.実際に確率的に示すことによって,得られた差異がどの程度珍しいのか,ということを示しています.例えば,p値が0.03,すなわち3%であれば,帰無仮説が正しいとした時に今得られた結果は3%でしか観察できないような珍しいことが起こっていることを示しています.こんなに珍しいことが起こったのは,その帰無仮説が正しくないからであり対立仮説を選ぼう!という論理のもとに対立仮説を採択することになります.
ここでは,t値とp値の計算方法については別書に譲ることとして,ざっくりとした理解で先に行きましょう.
t検定とは2群の「平均値」を比較する方法です.しかし,実はこれも一般線形モデルの枠組みの中で考えることが出来ます.ここではその考え方について説明します.そこには「ダミー変数」という考え方が必要になります.
一般線形モデルではこんなモデル式から考える,というような話をしたかと思います.
\[Y_i=\beta_1 X_1 + \alpha+\epsilon_i \]
数式の\(X_1\)に0を代入しましょう.
\[Y_i=\alpha+\epsilon_i \]
数式の\(X_1\)に1を代入しましょう.
\[Y_i=\beta_1 + \alpha+\epsilon_i \]
このように,0か1の数字を入れてあげると0を入れられたグループと1を割り振られたグループでどれだけ差があるのか,ということを評価することができます.
そして,その「差」がどの程度あるのかも比較することができます.ここでは,主観的幸福度に未婚者と既婚者の間に差があるのか否かを,先ほどと同じような流れで考えていきましょう.
t検定に当たるのは2つの群に差があるのか否か,です.「差がある」を対立仮説,「差があるとはいえない」を帰無仮説とします.したがって,以下のような仮説を立てることが出来ます.
はじめに,分析対象となるデータを読み込んでおきましょう.
以前紹介した“esquisse”を使っていただいて構いません.動画をご確認ください.
また,別の方法としてggplotguiを使ったプロットの方法についても紹介したいと思います.
library(ggplotgui) ggplot_shiny(exdataset)
そうすると新しいウィンドウが開きます.
以下の通りの作業をしましょう.
# You need the following package(s):
library("ggplot2")
# The code below will generate the graph:
graph <- ggplot(exdataset, aes(x = MAR, y = SUB_HAP)) +
geom_point(stat = 'summary', fun.y = 'mean') +
geom_errorbar(stat = 'summary', fun.data = 'mean_se',
width=0, fun.args = list(mult = 1.96)) +
theme_bw()
graph
# If you want the plot to be interactive,
# you need the following package(s):
library("plotly")
ggplotly(graph)
これも同様に,本当に差があるのかどうかは,感覚的には明らかになっても科学的な根拠がありません.同じように検定をして確かめてみましょう.
#"hapsat_model"というオブジェクトに,分析モデルを代入する. marhap_model<-lm(SUB_HAP ~ MAR, data = exdataset)
#分析結果の要約を出力する summary(marhap_model)
## ## Call: ## lm(formula = SUB_HAP ~ MAR, data = exdataset) ## ## Residuals: ## Min 1Q Median 3Q Max ## -6.6538 -1.6538 0.3462 1.3462 4.9391 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 6.65378 0.09274 71.74 <2e-16 *** ## MARNotMarried -1.59286 0.14499 -10.99 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 2.212 on 961 degrees of freedom ## Multiple R-squared: 0.1116, Adjusted R-squared: 0.1106 ## F-statistic: 120.7 on 1 and 961 DF, p-value: < 2.2e-16
## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 5.0609 0.1115 45.41 <2e-16 *** ## MAR 1.5929 0.1450 10.99 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
この結果について,またモデル式と共に説明します.この結果は\(\alpha\) が5.0609,\(\beta\) が1.5929という結果でした.したがって,モデル式は以下のように示すことができます.
\[Y_i=1.59291 X_1 + 5.0609+\epsilon_i \]
まずは係数について説明します.これは未婚者の場合と既婚者の場合について考えてみましょう.
未婚者の場合は\(X_1\)が0でした.したがって,以下のように示されます.
\[Y_i= 5.0609+\epsilon_i \]
既婚者の場合は\(X_1\)が1でした.したがって,以下のように示されます.
\[Y_i=1.59291 + 5.0609+\epsilon_i \]
library(huxtable) huxreg(marhap_model)
| (1) | |
|---|---|
| (Intercept) | 6.654 *** |
| (0.093) | |
| MARNotMarried | -1.593 *** |
| (0.145) | |
| N | 963 |
| R2 | 0.112 |
| logLik | -2130.084 |
| AIC | 4266.168 |
| *** p < 0.001; ** p < 0.01; * p < 0.05. | |
# 初めて使う時はinstall.packages("coefplot")が必要です.
library(coefplot)
coefplot(marhap_model)
library(stargazer) stargazer(marhap_model, type = "html", align=TRUE, title = "分析結果", out = "marhap_model.xls")
ダミー回帰分析モデルによって未婚者に比べて,既婚者の方が主観的幸福度が1.59高いこと0.001%水準で示された. (一緒に表を見せると良い.)
ダミー回帰分析モデルによって未婚者に比べて,既婚者の方が主観的幸福度が1.59高いことが示された.(t(961)=10.99, p<.001).
\((主観的幸福度)=1.59291 (t=10.99) ×(既婚ダミー) \\+ 5.0609 + \epsilon_i\)
今までは一般線形モデルの枠組みからt検定の紹介を,すなわちダミー回帰分析の1つとしてのt検定を紹介しました.一方で,普通のt検定は以下のように行うことができます.
最近はt検定にもいろいろな方法が提案されています.従来は等分散性を検定するF検定を実施し後に,等分散を仮定したスチューデント(Student)のt検定を行ったり,不等分散を仮定したウェルチ(Welch)のt検定を実施する,ということが行われてきました.
しかしながら,2回検定を行うことは「検定の多重性」の観点から問題ではないか,という指摘もあったりします.
そこで,最近ではF検定を実施せずに いきなりウェルチのt検定を行うことが多くなっています.
welch_t.testmodel<-t.test(SUB_HAP ~ MAR, data = exdataset) welch_t.testmodel
## ## Welch Two Sample t-test ## ## data: SUB_HAP by MAR ## t = 10.854, df = 808.29, p-value < 2.2e-16 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## 1.30479 1.88094 ## sample estimates: ## mean in group Married mean in group NotMarried ## 6.653779 5.060914
student_t.testmodel<-t.test(SUB_HAP ~ MAR, data = exdataset, var.equal = T) student_t.testmodel
## ## Two Sample t-test ## ## data: SUB_HAP by MAR ## t = 10.986, df = 961, p-value < 2.2e-16 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## 1.308323 1.877406 ## sample estimates: ## mean in group Married mean in group NotMarried ## 6.653779 5.060914
ちなみに,スチューデントのt検定と一般線形モデルにおけるダミー回帰モデルは結果が一致します.これは一般線形モデルが等分散性を仮定していることによります.
# モデルを作る
オブジェクト <- lm(応答変数 ~ 説明変数,
data = データセットの名前)
# 結果を出力する
summary(オブジェクト)
オブジェクト <- t.test(応答変数(数量データ) ~
説明変数(2つで分けられるもの),
data = データセットの名前)
オブジェクト
“CHI”は子どもの有無を尋ねる項目である.
これについて,以下の3つの分析を実施せよ.また,それぞれについてグラフも作成せよ.
* 子の有無による主観的幸福度の差を分析せよ.
* 子の有無による生活満足度の差を分析せよ.
* 子の有無による睡眠満足度の差を分析せよ.
分散分析とは,「3群以上の分散に差があるかどうか」を比較・分析するための方法です.その後「多重比較」という手法を用いて,「3群以上の平均値の差があるかどうか」を明らかにします.この授業では「1元配置分散分析」および「2元配置分散分析」というものについて説明します.いずれについても,説明変数が因子データ,応答変数が数値データとなります.
1元配置分散分析:「地域によって,主観的幸福度の分散・平均値が異なる」などのような,1つの要因によって影響を受けるかどうかを分析する手法です.
2元配置分散分析:「地域と未婚・既婚によって分散・平均値が主観的幸福度が異なる」,「地域と子の有無によって主観的幸福度が異なる」などのような,2つの要因によって影響を受けるかどうかを分析する手法です.
分散分析を一般線形モデルの枠組みで説明すると,平均値の推定がベースとなりますが,以下のように理解することができます.ここでは,「3つの群の影響を受ける」場合について,モデル式を元に説明します.また,以下では「分散分析モデル」という表現をします.
\[ Y_i=\beta_1 X_1 +\beta_2 X_2 + \alpha+\epsilon_i \]
このモデルでは\(X_1\)と \(X_2\) はそれぞれ(1, 0)の値を取る「ダミー変数」です.しかし,これでは\(\beta\) が2つしかありません.しかし,これだけで3つの群を表すことができます.以下には3つの条件についてモデル式を書き入れてあげたいと思います.
このモデルについて,平均値が異なるかどうかを調べます.特に,分散分析の場合は「分散分析表」と呼ばれるものを出して評価してあげます.
\[ Y_i=\beta_1 X_1 +\beta_2 X_2 + \alpha+\epsilon_i \]
\(X_1 =1\)と\(X_2 =0\):Bクラス \[ Y_i=\beta_1 +\alpha+\epsilon_i \]
\(X_1 =0\)と\(X_2 =1\):Cクラス \[ Y_i=\beta_2 + \alpha+\epsilon_i \]
\(X_1 =0\)と\(X_2 =0\):Aクラス \[ Y_i=\alpha+\epsilon_i \]
このモデル式からわかること:Aクラスに比べてBクラス/Cクラスの得点が高いか低いか
さて,それでは仮説を立ててみましょう.今回分析するテーマは「主観的幸福度(SUB_HAP)が地域(SUB_ARE)によって異なる」かどうかを分析します.一要因分散分析の場合は以下のような仮説を立てます.
この2つの仮説のもとに分析を行ないます.
今回の分析には,以下のモデルを前提とします. \[(主観的幸福度)=\beta_1 (北海道ダミー) +\beta_2 (東北ダミー) +\] \[\beta_3 (中部ダミー)+\beta_4 (近畿ダミー) +\beta_5 (中国ダミー)+\] \[\beta_6 (四国ダミー)+\beta_7 (九州ダミー) +\alpha+\epsilon_i \]
なお,このモデルではそれぞれの値は1か0の値しか取りません.
ex.東北地方のデータである場合には,東北ダミーが1,それ以外のダミー変数は0を取ります.
また,すべてのダミー変数が0の場合はコントロール群となる関東地方の値を示しています.
さて,例によってggplotguiを使いましょう.
以下のコードはConsole(コンソール)に直接打ち込みます.
library(ggplotgui) ggplot_shiny()
そうすると新しいウィンドウが開きます.
以下の通りの作業をしましょう.
# You need the following package(s):
library("ggplot2")
# The code below will generate the graph:
graph <- ggplot(exdataset, aes(x = ARE, y = SUB_HAP)) +
geom_point(stat = 'summary', fun.y = 'mean') +
geom_errorbar(stat = 'summary', fun.data = 'mean_se',
width=0, fun.args = list(mult = 1.96)) +
theme_bw()
そうすると,こんなグラフが算出されます.
graph
このグラフを見る限り,地域ごとに差があるかどうかはわかりません.以前,平均値を算出してみたことがありましたが,今回はそれぞれが「統計的に差がある」ということが言えるかどうかを考えたいと思います.
さて,分散分析モデルを作成してみましょう.
#"arehap_model"というオブジェクトに,分析モデルを代入する. arehap_model<-lm(SUB_HAP ~ ARE, data = exdataset)
#分析結果の要約を出力する summary(arehap_model)
## ## Call: ## lm(formula = SUB_HAP ~ ARE, data = exdataset) ## ## Residuals: ## Min 1Q Median 3Q Max ## -6.5429 -1.4308 0.1515 1.9043 4.7813 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 5.858108 0.192158 30.486 <2e-16 *** ## AREChugoku 0.572661 0.347849 1.646 0.1000 ## AREHokkaido 0.684749 0.439390 1.558 0.1195 ## AREKanto 0.237637 0.226845 1.048 0.2951 ## AREKinki -0.009623 0.264660 -0.036 0.9710 ## AREKyushu 0.228848 0.310363 0.737 0.4611 ## AREShikoku 0.530781 0.583547 0.910 0.3633 ## ARETohoku -0.639358 0.349733 -1.828 0.0678 . ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 2.338 on 955 degrees of freedom ## Multiple R-squared: 0.01418, Adjusted R-squared: 0.006954 ## F-statistic: 1.962 on 7 and 955 DF, p-value: 0.05729
Coefficientsだけ示します.Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.095745 0.120558 50.563 < 2e-16 ***
AREHokkaido 0.447112 0.413125 1.082 0.27941
ARETohoku -0.876995 0.316105 -2.774 0.00564 **
AREChubu -0.237637 0.226845 -1.048 0.29510
AREKinki -0.247260 0.218299 -1.133 0.25764
AREChugoku 0.335025 0.314020 1.067 0.28629
AREShikoku 0.293144 0.564036 0.520 0.60338
AREKyushu -0.008788 0.271909 -0.032 0.97422
6.095745である.\[ (主観的幸福度)=0.447112*(北海道)-0.876995*(東北)- \] \[ 0.237637*(中部)-0.247260*(近畿)+ \] \[ 0.335025*(中国)+0.293144*(四国)- \] \[ 0.008788*(九州)+6.095745 + \epsilon_i \]
#分散分析表 anova(arehap_model)
| Df | Sum Sq | Mean Sq | F value | Pr(>F) |
|---|---|---|---|---|
| 7 | 75.1 | 10.7 | 1.96 | 0.0573 |
| 955 | 5.22e+03 | 5.46 |
Analysis of Variance Table
Response: SUB_HAP
Df Sum Sq Mean Sq F value Pr(>F)
ARE 7 75.1 10.7238 1.9623 0.05729 .
Residuals 955 5218.9 5.4648
Dfは自由度を示しています.Sum Sqは平方和Mean Sqは平均平方F valueはF値Pr(>F)はp値を示しています.Response: SUB_HAP
Df Sum Sq Mean Sq F value Pr(>F)
ARE 7 75.1 10.7238 1.9623 0.05729 .
Residuals 955 5218.9 5.4648
SUB_HAPです.AREの自由度(分子自由度)は7:全部で8地域ある→N-1が自由度
DFはDegree of FreedomAREのF値は1.9623,P値は0.05729Residualsの自由度(分母自由度)は955:全部で963個のデータがあり,モデル式の\(\beta\)(パラメータ)で7つ,さらにもう1地域(=\(\alpha\)で使われる)を引いたもの.--- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
***で表される.**で表される.*で表される..で表される.library(huxtable) huxreg(arehap_model)
| (1) | |
|---|---|
| (Intercept) | 5.858 *** |
| (0.192) | |
| AREChugoku | 0.573 |
| (0.348) | |
| AREHokkaido | 0.685 |
| (0.439) | |
| AREKanto | 0.238 |
| (0.227) | |
| AREKinki | -0.010 |
| (0.265) | |
| AREKyushu | 0.229 |
| (0.310) | |
| AREShikoku | 0.531 |
| (0.584) | |
| ARETohoku | -0.639 |
| (0.350) | |
| N | 963 |
| R2 | 0.014 |
| logLik | -2180.170 |
| AIC | 4378.340 |
| *** p < 0.001; ** p < 0.01; * p < 0.05. | |
library(coefplot) coefplot(arehap_model)
library(stargazer) stargazer(arehap_model, type = "html", align=TRUE, title = "分析結果", out = "arehap_model.xls")
オブジェクト<-lm(応答変数 ~ 説明変数,
data = データセットの名前)
これについて,回帰分析/t検定の時は以下のコードを使っています.
summary(オブジェクト)
これについて,分散分析の時は以下のコードを使っています.
anova(オブジェクト)
“SUB_SAT”は生活満足度,“SUB_SLP”は睡眠満足度に関するデータであった(各10点尺度).これらを応答変数,地域を表す“ARE”を説明変数として,以下の2つの分析を実施せよ.
“SUB_SAT”は生活満足度,“SUB_SLP”は睡眠満足度に関するデータであった(各10点尺度).これらを応答変数,年代を表す“GEN”を説明変数として,以下の2つの分析を実施せよ.